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Abstract. The role of forest heterogeneity in the long-term, large-scale dynamics 
of forest fires is investigated by means of a cellular automata model and mean field 
approximation. Heterogeneity was conceived as trees (or acres of forest) with distinct 
strengths of resistance to burn. The scaling analysis of fire-size and fire-lifetime 
frequency distributions in the non-interacting fire steady-state limit indicates the 
breakdown of the power-law behavior whenever the resistance strength parameter 
R exceeds a certain value. For higher resistant strength, exponential behavior 
characterizes the frequency distributions, while power-law like behavior was observed 
for the lower resistant case in the same manner as reported in the literature for a 
homogeneous counterpart model. For the intermediate resistance strength, however, 
it may be described either by a stretched exponential or by a power-law plot whenever 
the fraction of recovering empty cells by susceptible trees not-exceeds or exceeds a 
certain threshold respectively, also suggesting a dynamical percolation transition with 
respect to the stationary forest density. 
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1. Introduction 

A forest-fire model was conceived by Bak, Chen and Tang (BCT)[T] in an attempt to 
study turbulence using the concepts of self- organized criticality (S0C)[2l |3]. In such 
a model, trees are supposed to grow on empty cells at a constant rate p while fire, 
once started by a single lightning strike, spreads through the system by deterministic 
local rules. The BCT model however, does not exhibit SOC, but rather a deterministic 
behavior marked by nontrivial spiral like fire fronts with a well defined length scale (1/p) 
within the limit p ^ [Ij. Drossel and Schwabl (DS) O |6] added a new ingredient to 
the BCT model with the aim of providing a mechanism to ignite small forest clusters in 
the same manner as the large ones. This new ingredient, the probability of a lightning 
strike /, led to power-law like behavior for both the fire-size (the total number of trees 
burned in a fire) and the fire-lifetime (the duration of fires in unities of time steps) 
frequency distributions, which is considered in the literature the signature of the SOC 
state within the limit f/p 0. However, such power-law distributions have been 
criticized in the literature[71 |8]. Nevertheless, the DS forest-fire model has become a 
paradigm for non- conservative systems exhibiting self-organized criticality [3]. 

By the end of the past century some authors began to perform statistics of burned 
areas of long series of events recorded in historical wildfire database with aim to find 
signatures of self-organized criticality (power-law) as firstly emphasized in references 
[9l [To] and later investigated in [HI [121 [13] • However, other frequency distributions of 
real fires were found in disagreement with a single power-law behavior being smooth 
but more approximating described by multiple power-law ranges delimited by cutoffs 
[m [m [151 [IS]. Recently, Corral and co-authors [17] found that the probability 
densities of forest-fire waiting times, defined as the time between successive fires above a 
certain size, verified a scaling law despite the fact that the associated fire-size frequency 
distribution does not follow a power-law, rather being fitted by a lognormal distribution. 

Concerning modeling, more recently the DS model has been explored and 
generalized by adding structural features and processes to the model in attempt to 
explain the fire-size multiple power-law or hump-shaped distributions observed in many 
forest-fire data [18]. Other distinct model was also proposed with the same aim by 
including the main characteristics vegetational growth and fire ignition and propagation 

m- 

In the present work, the DS forest-fire model was generalized by introducing 
heterogeneity into the forested environment, which consists of trees (or acres of forest) 
with more resistance to burn in comparison to the plain susceptible tree population. This 
heterogeneity within the tree populations - with regard to the strength of resistance to 
burn - should be more appropriated for describing fires in natural forests. Moreover, the 
present simulations also illustrate that such heterogeneity breaks the scaling behavior 
that appears in the homogeneous DS forest-fire model when the resistance strength to 
burn reaches a certain threshold. More specifically, the cellular automata (CA) DS 
forest-fire model [5] was generalized considering the simplest case of a forest composed 
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of two populations of trees - the resistant and the susceptible and which are 
randomly distributed on a two-dimensional square lattice. In addition, each CA- 
cell is considered to have a neighborhood comprising its eight surrounding cells so 
called Moore neighborhood in contrast with that composed by the first four nearest- 
neighbors cells (called von Neumann neighborhood) adopted by the DS model. Resistant 
trees are assumed to burn by contact whenever there are at least R burning trees 
in their neighborhood, while susceptible trees require at least a single neighboring 
burning tree in order to ignite in the next time step. Hence, R is in some sense 
the parameter that controls the strength of resistant to burn of the resistant-tree 
population. Such a deterministic mechanism differs from the stochastic immunity 
previously considered by Drossel and Schwabl [20], where the immunity to burn was 
defined by a certain probability controlling the ignition of a tree surrounded by burning 
trees. In the present model, the strength of resistance to burn has a deterministic 
character reinforcing the inherent local inhomogeneities of the distribution of trees along 
the forested environment. 

The coming sections are organized as follows: section [2] describes de cellular 
automata model and the limit conditions of the simulations, section [3] accounts for 
the study of the fire-size and time-size probability distributions focusing on their long 
tail behavior more precisely on the breaking of SOC, while section |4] is dedicated to 
study of the long term steady-state of the populations with special attention to the 
dynamical percolation transition. Section [5] summarizes the results and conclusions. 

2. Cellular automata model and simulations 
2.1. Cellular automata rules and parameters 

The present cellular automata model is defined on a L x L square lattice with periodic 
boundary conditions and Moore neighborhood. The variable associated to cellular 
automaton cells may assume one of the four possible states: empty cell (S), susceptible 
tree (S), resistant tree {TZ), or burning tree (B). Simulations start from a random 
configuration composed by susceptible and resistant trees, and empty cells. Afterwards, 
all cell states are synchronously up-dated according to the following rules: 

1. A tree grows on an empty cell with probability p, being susceptible with probability 
q or resistant with probability (1 — g). 

2. A susceptible or resistant tree is hit by a lightning strike with probability /, 
becoming a burning tree. 

3. A burning tree dies, becoming an empty cell. 

4. A susceptible tree with at least one burning neighboring tree becomes a burning 
tree. 

5. A resistant tree with at least R burning trees on its neighborhood becomes a burning 
tree. 

One time step is defined as the total time to update in parallel the entire lattice according 
to rules (l)-(5). Rules (1) and (2) play the role of the environment, which may be 
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associated with factors external to the intrinsic dynamic of fires such as hghtning and 
capacity for growth of trees, while rules (3) through (5) determine the microscopic 
dynamics of fire propagation. 

The parameter p is the probability that one tree growths on a given empty cell in one 
time step, while / is the probability that a lightning strikes a single tree {S or TV) in one 
time step. The latter represents external perturbations, which starts fires (avalanches) 
allowing that small and isolated forest clusters ignite sustaining fire at long-term scale. 
Both stochastic parameters p and / represent probabilities per time and grid cell and 
play the same role as the corresponding ones defined in the homogeneous DS model 

Concerning the tree population heterogeneity the new parameters introduced in the 
present model was the strength of resistance to burn R and the replenishment parameter 
q. According to the Moore neighborhood R may vary from 1 to 8. The trivial case R = 1 
characterizes the susceptible-tree S population (as in the DS model) and R chosen within 
the interval 2 to 8 characterizes the resistant-tree TZ population. Therefore the present 
CA model has two populations of trees: one composed of susceptible trees and the 
other composed with resistant trees all of them with the same strength of resistant to 
burn R. The parameter q fixed the fraction of new trees that grow on empty cells in 
the susceptible S state, hence the (1 — g) is the fraction of new trees that grow in the 
resistant TZ state. The long-term behavior of the forest population when R and q are 
varied is the focus of interest of this work. 

2.2. Time scales and power-law regime conditions 

The two dynamical characteristic time scales arise through the inverse of probabilities p 
and /. 1// corresponds to the average time between two consecutive lightning strikes, 
while 1/p corresponds to the average time between the birth of two consecutive trees. 
These two time scales are quite distinct in the limit f/p 0, with f ^ 0. At such a 
limit the forest has enough time to relax, since a huge amount of trees of order oi p/ f 
may grow during the time interval of 1// time steps before a new lightning strike. The 
limit p ^ also means that fires are extinguished before new trees are able to grow, thus 
preventing new trees from growing along the edges of a burning cluster before the entire 
cluster is burned down. On the other hand, the limit / — J- means that fires caused by 
different lightning strikes do not overlap. When considering such a limit regime the rules 
of the model can be adapted to reduce the simulation time enormously. Starting from 
an arbitrary initial configuration composed of empty, susceptible and resistant cells the 
successive configurations are generated by choosing a cell at random and 

1'. If the cell has a tree {S or TZ) the cell becomes a burning tree {B). Fire then spreads 

over the entire forest cluster to which it belongs, according to the previous rules 

(3)-(5), being parallel updated until the fire is extinguished. 
2'. If the cell is empty (£), p/ f new cells are chosen at random, assigning to the 

new empty cell one susceptible tree with probability q or one resistant tree with 
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probability (1 — q)- 

Notice that according to rule 1' the fire spreading occurs only within the cluster of 
trees hit by the lightning, hence the synchronous up-date of the whole lattice becomes 
restricted to the cluster, saving computation time. Moreover, rule 2' introduces a short 
cut in the time scale contracting the process of replenishment, also saving computation 
time. Such an alternative particular model was conceived to study the directed 
percolation model as considered in references [21] and [22]. In such a model version, 
apart from q and R, the relevant parameter is p/ f, the ratio of the growth and sparking 
rate [23], the limit described above is reached when p// — > oo. 

2.3. Studied cases 

The above described alternative model was adopted in the present paper to investigate 
the SOC regime f /p ^ with / 7^ 0. The simulations revealed that the system reaches a 
robust average steady-state configuration after a long transient regardless of the initial 
configuration. Distinct initial configurations were considered, including the extremes 
included in all cells at susceptible, resistant, or empty states. The final robust steady 
state was only reached after different transient times. 

Firstly the condition of equiprobable replenishment of trees by susceptible and 
resistant ones (parameter q = 0.5) was considered to study the particular cases were 
R =1, 2, 3 and 4. Under such conditions the fire-size and fire-lifetime frequency 
distributions were calculated after the system reached the steady state, which is 
characterized by constant average fractions of £, S and TZ cell states. 

The emerging results of the simulations reproduced the scaling power-law behavior 
presented by the fire-size frequency distribution found in the homogeneous DS model 
whenever the resistance strength R=l and 2. Such behavior found for R = 1 was 
expected since in this case both populations are susceptible actually corresponding to 
a homogeneous forest. However, the power- law behavior found for R = 2 may be 
interpreted as resulting from fire propagation in a heterogeneous forested environment 
but composed by susceptible and weak resistant trees. On the other hand, when the 
resistance strength was chosen R =4 and beyond, the overall behavior departed from 
such a scenario and the cumulative frequency distributions exhibits clear exponential 
decay, breaking, as expected, the presumed SOC like behavior. However, for the 
particular case when the strength of resistant of the TZ cell population was set equal 
-R = 3 the system presented a not expected intermediate behavior, the cumulative 
frequency distribution functions being fitted neither by a power-law or exponential 
curves. For instance, a stretched exponential function seems to be a good fitting, as 
presented and discussed in section [31 

Secondly, the variation of the replenishment parameter q within the interval [0, 1] 
was considered in the study of the particular cases of tree resistant population with 
resistant strength R =2, 3 or 4. Such study is detailed presented and discussed in section 
m It is worth to anticipate that for the intermediate case R = 3 the fire-size frequency 
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distribution tails can exhibit power-law decay whenever the fraction q of replenishment 
by S trees exceed a certain threshold Qc — 0.82. Such a change in the cumulative fire- 
size frequency distribution is followed by a change in the absorbing density steady state 
indicating the occurrence of a dynamical percolation phase-transition. 

3. Fire-size and fire-lifetime frequency distributions 

Simulations were carried out using the four-tap algorithm described by Ziff [21] and 
suggested by Grassberger [23] as the random number generator. The lattice sizes and 
number of fires however, were restricted by the CPU-time limitations, being L = 10, 000 
ioT R = l and i? = 2, and L = 20, 000 for i? > 3 cases, both with p/f = 10^ Periodic 
boundary conditions [torus conditions) were adopted. To save computer time, the 
initial configurations were chosen with 50% susceptible trees and 50% resistant chosen 
at random. Once the system had reached the steady-state, the size and lifetime of each 
of the fires were measured. The frequency distribution for fire-size p{s) (s being 
the size of the fires) and fire-lifetime p{t) [t being the fire duration) were estimated by 
the corresponding frequency histograms. Figure [1] illustrates the normalized frequency 
distribution of fire sizes for p/f = 10^ and R = 1 — 4. It can be observed that the 
plots for R = 1 and 2 display a power-law like behavior within almost three decades 
indicating that such a frequency distributions for the fire-size have long tails. 

In order to avoid strong fluctuations of the larger sizes caused by poor statistics, 
the normalized cumulative frequency distributions P{s) = p{s') ds'/ p{s') ds' and 
P{t) = j'^ pit') dt' I pit') dt' was also investigated. 

Figure |2] shows the plots of the cumulative fire-size frequency distribution P(s) for 
forests with a resistance strength oi R = 1—4, considering the fraction of replenishment 
of susceptible trees as g = 0.5. Dependence on such a parameter is discussed later on 
section m The plot for P{t) have similar characteristics and was skipped in the present 
manuscript. 

In the DS forest-fire model the relevant physical quantities are the fire-size and 
fire-lifetime frequency distributions. However, the critical behavior is expected when 
p// — )■ oo in the thermodynamic limit L — )■ oo. It is necessary, though, to take both 
limits, L oo, first, and p// — )■ oo afterwards, hence, the finite-size scaling analysis 
needs to be done in two variables, L and p/f. This feature was well studied in the 
literature. For instance, the scaling behavior shown by the DS forest-fire model reported 
in the literature (e.g. references [ZllHlEEj) is recovered in the present model when R = 1 
is established, in spite of considering the Moore neighborhood. For such a case the 
system represents a homogeneous forest, since only susceptible trees are considered, the 
q parameter becoming irrelevant. The plot for R = 2 displays the same qualitative 
behavior as that for R = 1, i.e. slowly decaying cumulative fire-size and fire-lifetime 
distributions are observed in accordance with (controversial) power-law P{s) ~ and 
P{t) ~ t^~^ behavior reported for the DS forest-fire model. Exponents r and (3 can be 
estimated by direct linear fitting. We found r = 1.18126(3) for i? = 1 and r = 1.18432(4) 
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Figure 1. (color on-line) Log-Log plot of the normalized fire-size frequency 
distributions for: weak resistant case i? = 2 (black), intermediate resistant case i? = 3 
(blue) and strong resistant case i? = 4 (green). Inset shows the homogeneous i? = 1 
case. Simulation parameters: p/f = 10^, q = 0.5. L = 10^, 10^ samples for i? = 1 
and 2 cases, and L = 2 x 10^, 10^ samples for i? = 3 and 4 cases. Many points were 
skipped in the R = I and 2 plots to avoid heavy figure file. 



for R = 2 when the scahng intervals are chosen by eyes, as shown in Figure [2j Analogous 
results may be obtained for the (3 exponent. Such exponents can be also estimated by the 
method of moments, introduced by De Menech et al [2S1 EZ], which has the advantage 
to consider the whole frequency distribution instead of choosing arbitrarily the cutoff 
by eyes. Such procedure was considered to study the homogeneous DS model in [25] . 
Adopting this procedure, the following scaling functions and associated exponents are 
defined for the cumulative fire-size and fire-lifetime frequency distributions, respectively. 



and 

where g{x) and h{x) are the corresponding scaling functions. The rate p/f has to go to 
infinity as L^^^. 
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Figure 2. Cumulative fire-size distributions for R ~ 1,2,3 and 4. Simulation 
parameters are: L = 10^ for i? = 1 and R = 2; L = 2 x 10^ for i? = 3 and i? = 4; 
p/ f = 10^, and q = 0.5. Dashed lines indicate the corresponding fittings. 

In the present study, the main motivation is to understand how heterogeneity can 
affect the large scale behavior of fires, so simulations were performed in lattices with 
fixed large enough sizes to accomplish this goal, doing the finite-size scaling analysis in 
the variable p/f in the range 10^ to 10^, similar to the method developed in reference 
[26|[27]. Concerning the cases where power-law behavior were observed, Tabled] displays 
the estimates of the scaling exponents r and A, (3 and defined by equation [1] for the 
present model with R = 1 and R = 2, as well as the corresponding exponents reported 
in reference [25] for homogeneous DS model, showing good agreement. 

The differences between both estimates may arise due to the fact that the scaling 
method considers the entire distribution, avoiding the choice of arbitrary cut-offs. We 
remark that independently of the considered method, both R = 1 and R = 2 cases 
display similar power-law behavior to that observed in the homogeneous DS model 
(with von Neumann neighborhood). 

When the resistance parameter has become fixed at i? > 4 (g = 0.5), it is unlikely 
that resistant trees will ignite by contact. The exponential distribution of frequencies 
emerges clearly when the graph is performed on mono-log scale. We also checked visually 
through snap-shootings of the steady state that the forest evolves to a configuration with 
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r 


1.098 


1.069 


1.08 ±0.01 


A 


1.096 


1.062 


1.09 ±0.01 




1.247 


1.228 


1.27 ±0.01 




0.579 


0.567 


0.59 ±0.01 



Table 1. Scaling exponents associated with the fire-size and fire hfetime frequency 
distribution in the steady state for homogeneous forests {R = 1) and inhoniogeneous 
forests with low-resistant trees [R = 2). Simulation parameters were q = 0.5, 
L = 10,000, and p/ f ranging from 10"* to 10^ 

most resistant trees forming compact groves, which prevent the occurrence of large fires. 
In practice, they become burning trees only when directly hit by a lightning strike. As 
the system evolves, the entire population of trees becomes resistant to fire and the fire- 
size decreases. An exponential decay in the cumulative fire-size frequency distribution 
following 

P{s) ~ exp (-s/so). So = 12.36 (3) 

was found. The cumulative fire-size frequency distributions P{s) for forests with 
resistance strength i? > 4 are almost all the same as the case R = 4 and were omitted 
for the present discussion. 

The intermediate case occurs when R = 3 (g = 0.5) where resistant trees are not 
so strong, but can still defeat fire. In such a case, neither power-law nor exponential 
behavior were observed. To understand the behavior of the case R = 3 and how it differs 
from the case R = 1, 2 (power law) or the R > 4 cases (exponential tails) we observed 
that in the former case groves of susceptible trees tends to a sparse fractal structure 
as g qc, as we noted in from snap-shootings of the steady state. When lightning 
hits a fractal cluster of susceptible trees, fire should spread in an anomalous manner. 
The fire front within a cluster marks the boundary of the diffusion of many random 
self-avoiding walker (SAW) starting from the same point on a given given realization. 
For clusters with fractal structures close to the one of the percolation cluster the fire 
front (SAW's) explore the labyrinthine geometry of the giant cluster leading to a slow 
relaxation regime. Therefore, the long fractal tails may increase the probability of fires 
reaching sites from different sides at once, increasing the probabilities of larger fires. 
Stretch exponentials have been recently used to correct describe the relaxation function 
for fractal time random walk process in one dimension [28]. In the present R = 3 case, 
the P{s) distribution may be fitted for instance with a stretched exponential function, 
i.e. 

P{s)=Aexp[-{s/sor] (4) 

with A = 2.75±0.01, sq = 1.08 and n = 0.36. The behavior of the fire lifetime frequency 
distributions are qualitatively identical. The k value obtained is in good agreement 
with the value 1/3, expected to relaxation processes that occur when the dilution of the 
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Time (x 10 steps) 

Figure 3. Time evolution of the forest density for four different values of R. 
Parameters are: p/ f = 10^, and q = 0.5. Lattice size L = 2 x IC*. 

support set reaches the percolation concentration of the underlying lattice. Such general 
behaviour was observed in the relaxation of the memory function in a random walks on 
the diluted hypercube, which was accurately parametrized by a stretched exponential 
over several orders of magnitude [29| [30] . 

4. Time evolution and steady-state regime 

Long-term behavior of forest density 

This section is devoted to investigate the role of the replenishment parameter q on the 
long-term behavior of the trees and empty cells population densities. The time evolution 
patterns of such population densities vary according to the strength of resistance to burn. 

Figure [3] displays the time-evolution patterns of the forest density (susceptible + 
resistant trees) for four different values of R and replenishing parameter q = 0.5. We 
note that the time scale is measured in units of 1//, and the densities are estimated 
just after the fires have been extinguished. Therefore, the plots exhibit curves with 
discontinuities everywhere, the drops representing the fires avalanches, whose cumulative 
fire-size frequency distributions were displayed in Figure El 
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When analyzed on a coarse-grained scale the scaling cases i? = 1, 2 exhibit highly 
fluctuating patterns, with stronger fluctuation in the R = 1 case than the R = 2. The 
transient period of resistant cases R> 3 however, exhibit few strong oscillations at the 
beginning of the process followed by a smooth increase leading to a high level of the 
steady state concentration. For the sake of comparison, the initial conflguration in all 
simulations was taken to be the only one with all empty cells. Therefore, the early 
stage of the dynamics is strongly non-linear, characterized by the continuous growth of 
the forest up to a critical density followed by a cascade of large flres. After this, forest 
replenishment resumes, and the process repeats itself until it achieves the corresponding 
dynamic equilibrium. However, in all plots, after the characteristic transient period, the 
forest density would evolve towards the corresponding steady state mean value. The 
scaling cases exhibit strong dispersion (not completely seen in the patterns of Figure [3]), 
while resistant cases display a smooth increase. It is worth noticing however, that the 
average steady state densities are independent of the initial conflguration in all studied 
cases and depend upon R and q. Nevertheless, the long transient periods depend upon 
such parameters, the size of the lattice and the initial conflgurations. 

4-2. Equilibrium steady-state population densities 

The above results shown in Figure [3] were obtained considering the replenishment 
parameter q = 0.5. The role of such a parameter on the steady-state forest composition 
was further investigated. The population density of each tree state and empty cells 
were calculated in the steady-state regime, varying the replenishment parameter q in 
the interval [0, 1] for the resistance strength R =2, 3 and 4. The plots of these population 
densities as a function of the replenishment q are illustrated in the flrst three panels of 
Figure H] for the heterogeneous forests, corresponding to three distinct scenarios, which 
can be discriminated according to the resistance strength parameter. In all scenarios, 
however, when g = the density of trees reaches a trivial absorbing state, in which all 
cells are occupied by resistant trees, while for g = 1 the steady state is composed by 
patches of susceptible trees occupying a fraction ^ 0.27 ± 0.02 of cells, the remanning 
being empties. Such latter state also corresponds to the steady-state density limit 
of the homogeneous forest model {R = 1 case) with Moore neighborhood. For the 
corresponding DS model with von Neumann neighborhood the equilibrium forest density 
reaches ~ 0.4084 ± 0.0001 [22]. 

For R = 2, an abrupt decay of the forest density is observed as soon as q departs 
from zero reaching an equilibrium value that diminishes towards the q = 1 limit one as 
illustrated in the panel (a) of Figure HI 

For R = 3 (intermediate case) the forest density remains close to 1 as g varies 
from zero to a critical value of qc ~ 0.82. Such a steady state is characterized by 
a high concentration of resistant trees leading to large patches, which prevents the 
propagation of flre. Within this interval, the density of resistant (susceptible) trees 
decreases (increases) monotonically with q, while the density of the empty cells remains 
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Figure 4. Steady-state density of trees and empty cells as a function of the 
replenishing parameter g, for simulations with parameters p// = 10^, L = 10^ averaged 
over 10'^ events (error bars are smaller than the symbols as illustrated). Panels (a) 
i? = 2, (b) i? = 3, (c) i? = 4 , and (d) SOC limit in mean-field approximation. Insets 
in panels (b), (c) and (d) illustrate the amplification of the critical regions. Symbol 
labels: A resistant trees, Q susceptible, V empty cells , and □ forest density. Open 
symbols for < q < 1 and solid ones for the exact q = I limit case. In all cases 
the steady-state was reached by an automatic procedure based on the mean slope, 
averaged over 10^ simulations. 
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quite close to zero. However, when q > 0.82 an abrupt decay of the forest density 
occurs concomitantly with an increase of empty cells and susceptible trees. In the 
interval (0.82, 1] the forest density decays towards the homogeneous g = 1 (or i? = 1) 
limit value 0.27 ± 0.02. In this interval the density of resistant trees decreases to zero 
while the density of susceptible trees reaches a maximum, followed by a decay to the 
limit of 0.27 ± 0.02. Since for q = 1 only susceptible trees survive in the steady-state. 

For -R > 4, the steady-state equilibrium forest density reaches values almost close 
to one within the whole interval of values < g < 1, but abruptly decaying towards the 
limit of 0.27 ± 0.02 as g — 1, as illustrated in Figure H] panel (c). 

4-3. Mean- field approach 

The mean field approximation (MFA) has been applied to investigate the dynamics and 
steady state of propagation models, such as the forest fire [311 [321 133] and epidemiological 
models[33l[3ll|35] in complex networks. In this framework, the steady-state equilibrium 
forest density of the cellular- automata model defined by the above rules (l)-(5) is studied 
following the approach defined in reference [36] for a similar CA model. 

The fractions of each cell state Dr^ = (x = S, TZ, £ or B) can be obtained from 
the coupled difference equations for AD^ = D^it + 1) — D^it), i.e. 



ADe = DB{t)-pDe{t) 
ADs = qpDeit) - [6 + /(I - 6)] Dsit) 
ADn = (1 - q)pDe{t) - [^r + /(I - $/?,)] Dnit) 
ADb = [QDs{t) + ^RDn{t)] + f[{l-Q)Ds{t) + 
+ Dnit)] -Dsit) 



where 9(t) = 1 — {l — Dsit))^ is the probability that a S cell has at least one surrounding 
B cell in the time step t and 



is the probability that a TZ cell has at least R neighboring B cells at time step t. Such 
densities of cells must obey the normalization condition at each time step. 



In the steady-state D^ D^, that is, AD^ = ADg = ADn = ADq = 0, hence 
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'^nit) = E in) Dsity^l - Dsit)) 



8-n 



n=R 



Dsit) + Dnit) + Dsit) + Dsit) = 1 




Ds 



qDs 



(i-/)e + / 

il-q) Ds 



(1 -/)$« + / 
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where -D^ is given by the solutions of 



1 + - + 
P 



+ 



:i-/)e + / (1-/)$^ + / 

with G and $r being their corresponding values at the steady state. 

A similar plot can be drawn for the SOC limit condition studied here, say p — )■ 
/ ^ with f/p^O. For < g < 1 the fractions of burning trees Db, empty cells 
resistant D-ji and susceptible Dg trees are given by 



Ds 



X{q) f 



^0, 



(l-q) p 

7qX{q) 
56X{q) + 16(1 



f) 



1 + 32 (1 - g)/49] 



(5) 

(6) 
(7) 
(8) 

(9) 



For g = 1 however, the MFA-limit densities are 0, 7/8 and 1/8 for the resistant, empties 
and susceptible trees respectively, in slight numerical disagreement with the results 
obtained by simulation. The MFA results are displayed in the panel (d) of figure HI It 
is worth to mention that MFA-SOC limit is independent of the resistance strength, i?, 
as expected from the absence of local correlations in the MFA. 



5. Concluding remarks 

The present work investigates the role of heterogeneity in the previous homogeneous 
forest-fire model reported in the literature. Heterogeneity was considered by introducing 
a population of trees with a certain strength of resistant to burn defined by a local 
interaction and by a certain rate of replenishment. The main result emerging from the 
present study suggests that such ingredients play an essential role in the propagation 
process of forest-fires, changing the long-term, large-scale behavior of the overall 
phenomenon. The combination of both ingredients may lead the systems to reach 
two scenarios: a steady state characterized by power-law (long tail) fire-size and fire- 
lifetime frequency distributions, which is identified in the literature with the so-called 
self- organized critical state, or to converge to an exponential decay of such distribution, 
the latter leading to a forested environment composed by a background of resistant 
trees with patches of weak-resistant susceptible trees. For certain critical values of 
the parameters governing the heterogeneity the system steady state may undergoes a 
dynamical phase transition between absorbent states characterizing such scenarios. 

Real forests are very complex non-equilibrium, non-linear systems and many 
factors contribute to heterogeneities across the landscape of forested land. A forest 
in equilibrium can be characterized by having a species composition with age structure 
unchanged for long periods of time. Probably, few temperate forests may be found near 
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an equilibrium state, since fire, recurring pests, human interventions and other factors all 
tend to drive the system to a stochastic equilibrium. Data from real forest fires have been 
collected in many countries over the past decades and analyzed under many aspects. 
Regarding the frequency distribution of size of fires, many authors observed a clear 
power law behavior, while others have found deviations from this behavior but without 
doing a careful analysis of the diversity, age structure and other special characteristics of 
tree populations in the forests studied. The results of this study suggest that deviations 
in the frequency distribution of fire-sizes observed in many forests could be due to 
inhomogeneities in the populations of trees with respect to their resistance strength to 
ignition and its rate of replenishment. Further, such distinct populations evolve to a 
critical state of stochastic equilibrium. Further studies to characterize the nature of this 
transition is now being considered. 
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